NASA Contractor Report 145357 


A SIMULATION MODEL OF THE 
PLANETARY BOUNDARY LAYER 
AT KENNEDY SPACE CENTER 

Research Laboratories, Inc. 

FVinceton, N J 08540 

B.C. Hwang 


Contract NAS1-14987 
March 1978 


IW\SA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 



1. Report No. NASA Contractor 2. Government Accession No. 

Report 145357 

3. Recipient's Catalog No. 

4. Title and Subtitle 

A SIMULATION MODEL OF THE PLANETARY BOUNDARY 
LAYER AT KENNEDY SPACE CENTER 

5. Report Date 

March 1978 

6. Performing Organization Code 

7. Author(s) 

BaoChuan Hwang 

8. Performing Organization Report No. 

AeroChem TP-369 

10. Work Unit No. 

9. Performing Organization Name and Address 

AeroChem Research Laboratories, Inc. 
P.O. Box 12 
Princeton, NJ 08540 

11. Contract or Grant No. 

NASl-14987 

13. Type of Report and Period Covered 

Final - 24: June 1977 to 
28 March 1978 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 

14. Sponsoring Agency Code 


15, Supplementary Notes 


16. Abstract 


A simulation model which predicts the behavior of the Atmospheric 
Boundary Layer has been developed and coded. The model is partially 
evaluated by comparing it with laboratory measurements and the sounding 
measurements at Kennedy Space Center. The applicability of such an 
approach should prove quite widespread. 


18. Distribution Statement 

Unclassified - Unlimited 


page) 21. No. of Pages 22. Price* 

22 


17. Key Words (Suggested by Author(s)) 

Atmospheric Boundary Layer Simulation 

Second-Order Closure 

Turbulence 

20. Security Classif. (of this 
Unclassified 


19. Security Classif. (of this report) 

Unclassified 


For sale by the National Technical Information Service, Springfield, Virginia 22161 























TABLE OF CONTENTS 

Page 

SUMMARY ■ 1 

I. INTRODUCTION 2 

II. SIMULATED ATMOSPHERIC BOUNDARY LAYER 3 

A. The Basic Equations 3 

B. Modeling Assumptions 5 

C. One-Dimensional Model 7 

D. The Empirical Constants and Length Scales 8 

E. Initial and Boundary Conditions 9 

F. The Grid System and Finite Difference Approximation 11 

III. MODEL PREDICTIONS 15 

A. Laboratory Data Comparisons 15 

B. Sounding Measurement Comparisons 15 

IV. DISCUSSION AND CONCLUSIONS 18 

V. REFERENCES 20 


ill 



SUMMARY 


A model simulating the atmospheric boundary layer (ABL) at KSC is developed 
by using a simplified second-order closure scheme. This model predicts the time 
variations in the ABL of (i) wind velocity and temperature, (ii) inversion layer 
height, and (iii) turbulent Reynolds stress and temperature flux. Partial veri- 
fications of the model were made by comparing it with the wind tunnel measure- 
ments conducted by Arya and the regular sounding measurements at KSC on December 
10, 1974. Favorable correlations of the temperature profiles can be seen in 
both comparisons. A disagreement between the predicted wind velocity after 2.5 
hours and the actual sounding measurement at that time, in which the eastern 
component of the wind velocity is stronger than that measured, is caused by 
using a mandatory stronger stress near the ground as input. 

At present, the simulation model ignores (1) the horizontal advection and 
diffusion, (2) the buoyant force effect of the water within the ABL, and (3) 
the radiation effects. These neglected dynamic features may be Important in 
some meteorological conditions at KSC. Further investigation on the effects of 
these dynamic features is needed. 

The computer code of this model is stored in the computer facility at 
NASA/Langley. 
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I. INTRODUCTION 


The objective of this program is to simulate the atmospheric boundary 
layer (ABL) at Kennedy Space Center using a second-order closure scheme (or 
mean Rejmolds stress scheme) . Such simulation will provide the time variations 
in the ABL of (i) mean physical quantities such as wind velocity and temper- 
ature, (ii) inversion layer height, and (iii) mean turbulent Reynolds stress 
and temperature flux (characteristics of turbulence) . These physical /mete- 
orological quantities are required for: 

1. Assessing the environmental impact of a rocket launch several hours 
or even one or two days in advance. 

2. Input to diffusion models, e.g., NASA/MDM, for real-time calculations. 

3. Calibrating, locating, and orienting measurement instruments for 
remote sensing/measurements. 

4. Studying possible weather modifications due to rocket exhaust. 

5. Modeling the atmospheric chemistry in the ABL. 

6. Vehicle design and safety analysis. 

In addition, this Simulated Atmospheric Boundary Layer (SABL) model can act 
as a nesting model incorporated in the large-scale forecasting model (e.g., 
the Limited Area Time Mash Model in the National Weather Service) to forecast 
local weather. 

Progress of this program has been described previously in the first ^ and 
second^ quarterly progress reports. This report presents a summary of this 
work and discusses additional results obtained in the past two months. 

The support from Scott Wagner, Joseph Mathis, and Richard Bendura of 
NASA/Langley (Environmental Field Measurements Branch; Marine and Applications 
Technology Division) is gratefully acknowledged. 
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II . SIMULATED ATMOSPHERIC BOUNDARY LAYER 


A. THE BASIC EQUATIONS 


The basic equations employed in the simulation model use the 
following quantities; 


virtual temperature 


potential temperature 


T = (1 + 0.61 S)T 

V 




O . 2 SS 


and virtual potential temperature 


0 = 0(1 + 0.61 S) 

V 


where T is the absolute temperature, S the specific humidity, and P© the ref- 
erence pressure (1000 mb). In the above, the tilde refers to Instantaneous 
values, whereas below, the tilde will be removed to indicate mean quantities. 
Within the ABL, it is adequate to assume that the velocity change, due to 
variations in the density produced by variations of pressure (but not of 
temperature) ,_^is small. We also assume a constant (both time and spaclal 
independent) T© corresponding to 0©^, such that the density of air is 
described by 

*r p s Po + i 5** V©i s Po + P ©1 
\ 30 ^/ V V 


where 0i^ = ~ ©ov P = (^) the coefficients cf terminal expansion. 

The latter assumption is a strong Boussinesq approximation, permissible for 
dry air; it is a first-order approximation for including water vapor. In- 
corporating the above assumptions, the mean equations of motion for the 
velocity, U, V, and the virtual potential temperature, 0^, are 



8U 


3t 3xj^ ^jklVl 


• 5 ^ - g,p0 + vV*U. 

9Xj V j 


30 , 


5T + <\®v* 


( 1 ) 


( 2 ) 

(3) 
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where P is the mean kinematic pressure, gj = (o, o, — g) the gravitation vector, 
fj = (o, o, f) the Coriolis parameter, V the kinematic viscosity, and a the 
thermal diffusivity. The overbars represent ensemble averages; the lower case 
terms are the fluctuation components and are governed by (subtracting the equa- 
tion of motion from the mean Eqs. (1) to (3)) 


3x, 


3u 


3t ^j“k ^jkl^k“l' 

36 a 




Uj^0) = aV*0 


(4) 

(5) 

( 6 ) 


Since the mean Eqs. (1) to (3) involve the Reynolds stress u^Uj and the heat 

flux term uiQ, the mean equations are closed by Introducing the equations for 
the Reynolds stress and heat flux from Eqs. (5) and (6), 


3u.u. 

3t 
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U, 0^ + u, 0^ — a 
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3x, 
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B. MODELING ASSUMPTIONS 


Since the higher order terms. 


the diffusion terms u.u.u, and u.u,6 

1 3 k i j 


and the terms including the fluctuation of pressure, e.g., 

/ 3u. 9u . \ 

p Six^j called the energy redistribution terms), are involved 

'‘1 ^^^i' 


in Eqs. (7) to (9), these unknowns must be modeled if one does not intend to 
invoke a higher order closure scheme. The models adopted for these terms 
are based on the works of Rotta,® Kolmogoroff and Mellor and Herring.® They 
are 


/ 3u. 3u. \ 

1 i 

1 

= -_3_ 

3S.1 

( u.u. 

V 1 1 - 
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3x1. 
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— qXj 
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(17) 


u,e* 
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— qXs 


36 = 


* 

p6 = 0 


where q“ is twice the turbulent kinetic energy, u^ + + w" 

(10) to (18) into Eqs. (7) to (9) gives 


(18) 

Substituting Eqs. 
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8U _ 

~ ®"k ^ 
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D6^ 

Dt 
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k 3Xj^ A 2 


( 21 ) 


Although Eqs. (1) to (3) and (19) to (21) represent a closed second- 
order system, they are still too complicated to solve. The present program 
adopts a simplified model based on the so-called "level 3 model" by Mellor and 
Yamada, in which the off-diagonal turbulent diffusion and advection of Reynolds 
stress and flux terms were Ignored using the order of magnitude argument. Such 
a simplified model has been shown to give a favorable comparison with the 
Wangara experimental results.^ 
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C. QNE-DIMEN$1QNAL MODEL 


The one-dimensional simplified second-order closure system we plan 

to solve is 



( 22 ) 

( 23 ) 

( 24 ) 

( 25 ) 

( 26 ) 
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( 29 ) 


— — SXJ — — 3V — — 3y 

In the above expressions, = - uw Pyy = - vw — , and P^^ = - uw 

The Coriolis terms have been excluded in the turbulence equations; they may be 
shown to be small in the ABL. 

In practical geophysical fluid dynamics problems, the eddy viscosity 
and eddy conductivity are often used and are defined as 


t — — N X. 9V\ 

< - uw, - vw ) = % . 3^; 


(30) 


90 


w0 = Ki 


•H 3z 

Using Eqs. (27) to (29), and can be expressed as 

Km = 9 ^ % 




= q il S 


H 


where 


(31) 


_ 3Ai(w^ - Cxq=^)/q^ + 9AiAg£Pgw6/q^ 

~ 1 + 9AxAa)l2gg O0^/azy/q2 


,, ^ 2AgwVq^ 

H 1 + lAaBail^eg O0^/9z)q2 

and w^ is given by Eq. (27). 


D. THE EMPIRICAL CONSTANTS AND LENGTH SCALES 

The length scales (S-i, Ai, %z^ Az) incorporated in the model are 
assumed to be proportional to a master length scale I, i.e., 

(£i, Ai, Zz, Az) = (Ai , Bi, Az, Bz)!>' 


'8 



Although this assumption has not yet been verified, several studies by Mellor 
and Yamada® indicate that it may be viable. The constants A», Bi, A 2 , B 2 , 
and C should be determined from data. The present program, employing the 
values given by Mellor and Yamada,^ gives 

(Ai, Bi, A 2 , B 2 , C) = (0.74, 16.6, 0.92, 10.1, 0.08) 

These values were determined® by matching the turbulent characteristics of the 
measured neutral stratified wall data to the corresponding values in the model. 
Local equilibrium was assumed in interpreting the values obtained. 

The master scale, H, is determined by solving the equation 



(32) 

While one cannot place great confidence in Eq. (32), an alterna- 
tive for determining £ can be used»^ 


£ 


kz 



2.7 X 10““ V 

. j ^ 

where Vg is the geostrophic wind speed and k is the Von Karman constant. 

E. INITIAL AND BOUNDARY CONDITIONS. 

The initial value for the mean quantities U, V, and Oy for the one- 
dimensional model can be taken from direct rawinsonde data. Initial values 

for the turbulence moments q^ and 0* are given by invoking local equilibrium, 
i.e., that the total production balances the dissipation. Therefore we have, 
from Eq. (25) , 


and 


aT ■ “ If “ ™ 


(33) 




(34) 
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Solving the algebraic Eqs. (33), (34), (27), (28), (29) and again employing 
expressions (31) and (32) we obtain 


- (UW. V») . iqS„ 


( 35 ) 


where 


— w9 


30 


M 


3Ai 


Yi - C - (6Ai 3A2)r/Bi . _ p. 

‘ Yi - Yar + 3A,r/Bj ^ 


H 


= 3Aa(Yi - Yar) 


Yi = I - ( 2 Ai/Bi) 

Ya = (B2/B1) + (6A1/B1) 

r = Rj/(1 - Rj) 


(36) 


In these expressions, Is the flux Richardson number, l.e., 

Rf = - Pg^/ (?xx Pyy) 


The flux Richardson number can be related to gradient Richardson, R^ ■ 


^ (i)i 


as 


R. 








where P 


Sm 


’^t 


is the turbulent Prandtl number. 


relations 


The surface boundary conditions for U and V satisfy the neutral flow 


u 

(U, V) = (cos Cl sin C) (37) 

where u^ = [( — uw)^ + ( — w)^]^p is the surface friction velocity and C is 
the angle of the surface friction velocity. 

Based on Eq. (37) the following relation is used in the program for 
the surface condition of U and V 


Lio 



(Ux, Vi) 


= (Ua, Va) 


ln(zi/zo) 

InCza/zo) 


The surface boundary conditions for the 0y are set to be either 


©V “ “ 0) 


or 


K, 


H 



z=o 


(w0) 


z=0 


The surface boundary conditions for q* and 0^ determined from Ref. 8 are 
employed , i . e . , 


q^(0) 
^(0) = 


n 2/ 3 „ 2 

u. 


2//^n = (w9)^ . Ba , 


The upper boundary conditions for U, V, 0„, q, and 0* are 




0 



F 


and 

q*(h) = 02(h) = 0 

where F (K°m~^) is the observation value of the virtual potential temperature 
gradient . 


F. THE GRID SYSTEM AND FINITE DIFFERENCE APPROXIMATION 

Since (1) the master length scale H is proportional to z near the 
surface but tends toward a constant away from the surface and (2) an accurate 

Az 

computation requires a coordinate system such that -j- = const. , we choose a 
grid system {zjj} so that 

Az = Ain Zjj for lower portion of the ABL 

Az = AZjj for upper portion of the ABL 

Table I shows the 40 point distribution as an example of the chosen grid 
system. 
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TABLE I 


THE SELECTED FORTY POINT GRID SYSTEM 
Boundary Layer Height Assumed = 2500 m 


Grid Index 
n 

Height 
Z (m) 

Grid Index 
n 

Height 
Z (m) 

1 

0.0 

21 

1015.6 

2 

0.6 

22 

1093.8 

3 

1.2 

23 

1171.9 

4 

2.4 

24 

1250.0 

5 

4.9 

25 

1328.1 

6 

9.8 

26 

1406.3 

7 

19.5 

27 

1484.4 

8 

39.0 

28 

1562.5 

9 

78.1 

29 

1640.6 

10 

156.3 

30 

1718.8 

11 

234.4 

31 

1796.9 

12 

312.5 

32 

1875.0 

13 

390.6 

33 

1953.1 

14 

468.8 

34 

2031.3 

15 

546.9 

35 

2109.4 

16 

625.0 

36 

2187.5 

17 

703.1 

37 

2265.6 

18 

781.3 

38 

2343.8 

19 

859.4 

39 

2421.9 

20 

937.5 

40 

2500.0 
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The differential equations we plan to solve are (22) to (26) and 
(32) . They can be reduced to a general expression given by 


3 ({> 

3t 



- Pa If - PsCj) + 


Table II indicates (j>. Pi, P 2 , P 3 , 


P 4 for Eqs. (22) to (26) and (32). 


(38) 


TABLE II 

IDENTIFICATION OF THE VARIABLE <J) 
AND THE PARAMETERS Pi-Pfc IN EQ. (38) 


( 22 , 

23) 

(24) 


0 


Pi 




H 


I 2 . 

0 

W 


if 


0 


fT 


IfV, 


i'' ^ “ If h "r 


(25) 


a! 5 


qAi 


2#- 

Ai 


— 3U — 9V — 

uw ^ - vw — + pgw 0 


(26) ^ qX 3 


A» 


— w 6 


30 

V 

3z 


(32) 






•8 ^ fxx ^y + P 


J 


- ^ |l + 1.33 
B. I 



2 
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The leapfrog finite difference expression for Eq. (38) Is 


. k+1 . k-1 

2At 


+ 2 


i T J 




k 

j+1 


2(Az)2 ) J+1 

^ ♦i-1 [h j + 


1 

2Az 




k .k+1 
J 



(39) 


The finite difference expression (39) can be rewritten as 


BlfAl Oy • • • aO 


A *t+l 


+ 2At(P*)i'^" 

^2 f ^2 f ^2 t 0 • • *0 

• 


• 


• • 

• 


• 

s 

• • 

• 

• 

• 


• 


• • 

* S’ s 


. k+1 


+ 21t(P.)|; 


(AO) 


The matrix form (AO) can easily be solved by using the Thomas 
algorithm to give 


<> 


k+1 

N 


V 


N 


and 


4> 


k+1 


V4 + S. 


k+1 


J ^j+1 


for j < N 


where 


V» 



2At 



and 


'J 


- (Bi + C. S, ,) 


for 1 < J < N 


also 


S» 


- Ai 


and 



+ C. 



1 ^ j < N 
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III. MODEL PREDICTIONS 


The model has been partially evaluated by comparing it with two sets of 
data, the laboratory measurements carried out by Arya^° and the regular sound- 
ing measurements at KSC on December 10, 1974. They are discussed individually 
below. 

A. LABORATORY DATA COMPARISONS 


A wind tunnel simulation of a stably stratified surface layer was used 
for comparison because (1) it provides a complete data set including turbulence 
information and (2) the turbulence behavior of a laboratory simulated surface 
layer is less complex than that of an atmospheric boundary layer in nature. 

The temperature profile calculated using the model will be compared with the 
measured profile. Based on Eqs. (35) and (36), the temperature profile is 
related to the velocity profile by 


30 wQ 3U 

— = ' Pr — 

3z t 3z 


(41) 


Since the heat flux and the shear stress are approximately constant in the 
surface layer, the temperature profile can be easily obtained by integrating 
Eq. (41) with the given velocity profile and the turbulent Prandtl number. 

The turbulent Prandtl number in Eq. (41) is dependent upon the gradient Richard- 
son number as deduced from the model described in Section II. E. Using the 
turbulent Prandtl number from the model and the measured values for the veloc- 
ity profile, heat flux, and shear stress, a comparison was made of the calcu- 
lated and measured temperature profiles as shown in Fig. 1. The correlation 
appears very good. The comparison is restricted to between 0.2 and 0.8 Z^, 

where Z^ is the surface layer height, because the laboratory data show approxi- 
mately constant flux and constant stress in this range only. 

B. SOUNDING MEASUREMENT COMPARISONS 


Since the stationary state, the equilibrium state, the constant stress, 
arid the constant flux in the laboratory-simulated surface layer rarely exist 
in the natural ABL, the model is used to make a calculation and to compare it 
with the existing measured data in an ABL. Such comparisons will also guide 
further refinements of the model. A theoretical calculation using the model 
for the physical quantities in the ABL of KSC over a three hour period during 
the night of December 10, 1974 was made. The resultant wind velocities and 
temperatures are compared with the available sounding data in Fig. 2. 


A satisfactory agreement for the temperature and the southern component 
of the wind velocity, V, can be seen in the figure. On the other hand, the 
predicted value of the eastern component of the velocity, U, exceeds the 
measured value. This causes an error of about 20° in the predicted wind 
direction toward the east and a factor of 2 greater predicted than measured 
wind speed. Further discussion of this disagreement is given in the following 
section. 
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FIGURE 1 COMPARISON OF THE TEMPERATURE PROFILE BETWEEN 
THE WIND TUNNEL MEASUREMENTS*® AND THE SABL MODEL 
CALCULATIONS 


Kev 


Model calculation 


A Arya 




HEIGHT, 



TEMPERATURE, ®C VELOCITY, m sec"' 


FIGURE 2 COMPARISONS BETWEEN THE PREDICTED WIND VELOCITY AND TEMPERATURE AND THE 
ACTUAL SOUNDING MEASUREMENTS AT KSC ON 0441 Z, DECEMBER 10, 1974 

The sounding measurement on 0211 Z, December 10 is used in the calculation as input. 

Key: Predicted value; Actual measurement 





IV. DISCUSSION AND CONCLUSIONS 


At present, the simulation model which predicts the wind velocity and 
temperature throughout a horizontal quasi-homogeneous ABL has been developed, 
coded, and partially evaluated. This model is based on several simplifica- 
tions, viz., that (1) horizontal advection and horizontal diffusion are not 
important, (2) the dynamic effects of water (hvimidity) within the ABL on the 
buoyant mechanism can be neglected, and (3) radiation and the interactions of 
large phenomena, e.g., long wave propagations, do not significantly affect the 
predictions. Although these simplifications are not valid in some meteorolog- 
ical conditions, it is not now clear how much error each simplification can 
cause, and how important each one is under various meteorological patterns in 
the KSC area. In other words, appealing to Occam's razor, an extensive veri- 
fication of the developed model against various meteorological conditions at 
KSC should be made before attempting to further modify the code to include 
some of the ignored dynamic features. 

As is common in running a theoretical model for realistic problems, some 
difficulties exist in preparing the input information required to make pre- 
dictions. In this specific case, problem areas are the geostrophic wind at 
the surface, the ground temperature (or flux) history over the period covered 
by the prediction, and the stress near the ground. Although some established 
theorems to estimate these quantities do exist, these theorems unfortunately 
require more extensively measured data than are available in most situations. 
For example, the stress near the ground can be interpolated^^ from the local 
velocity profile near the ground, which is measured from the two tallest tow- 
ers in the KSC area. Since these two towers are located near the coastline 
and far from the sounding measurement site, it is doubtful that the stress 
interpolated from them is the same as that at the sounding measurement site. 

A higher value for this stress at the towers than at the sounding site can 
especially be expected in a land/ sea breeze condition because such a condition 
will intensify the turbulence strength and consequently the stress. As in the 
previous sounding comparisons, the input value of the stress near the ground 
was that interpolated from the tower measurements at as much as five hours 
after the time of the start of the calculations. The later (midnight) tower 
measurements were used because of the lack of tower data at the time desired. 
Since the land breeze is strongest at midnight, the stress used for the input 
is certainly higher than it should be; the eastern component of the stress is 
especially large because the sea is to the east at KSC. Therefore the eastern 
component of the predicted wind velocity in the comparison shown in Fig. 2 was 
stronger than the measured values because the higher mandatory input stress 
enhanced the vertical mixing and smoothed out the velocity differences between 
the higher and lower levels. Thus it is likely that the predicted value will 
compare better with the measurement data if more accurate surface stress is 
used as input. 

We had proposed earlier to evaluate the model by comparing it directly 
with the tower measurements at KSC and the Savannah River plant. Because 
previous studies® have shown that the Monin-Obukhov theory correlates well with 
tower data at KSC and the present model provides excellent agreement with the 
Monin-Obukhov predictions, the comparison between the model results and the 
KSC tower measurements need not be pursued. 
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In the present model, the initial value of the turbulent kinetic energy 
is obtained with the assumption that the total energy production (buoyant and 
shear production) balances the dissipation. This conflicts with the conclusion 
drawn from the spectra study of the horizontal turbulence at KSC by Fichtl and 
McVehil.^^ They have shown that the energy budget indicates a balance of 
dissipation and shear production only. Buoyant energy production apparently 
cancels the vertical flux divergence. Some uncertainties remaining in the work 
of Fichtl and McVehil should be checked before modifying the present model. 

For instance, the buoyant subrange was not clearly shown in their crude measured 
spectra. 

In conclusion, the present model simulating the ABL at Kennedy Space Center 
using a simplified second-order closure scheme demonstrates the potential to 
solve the very complex problem of predicting the behavior of the ABL. The 
applicability of such an approach should prove quite widespread; further veri- 
fication and/or refinement of the present model should be encouraged. 

The computer code for the simulation model is stored in the computer 
facility at NASA/Langley. 
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